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ABSTRACT 


An iteration method for solving the implicit difference equations 
associated with three-dimensional nonlinear parabolic differential equations 
is derived and analyzed. The method is applied to the high Reynolds number 
laminar viscous flow around a cone at high angle of attack, The requirements 

- which must be met to ensure convergence of the iterations are obtained, in 
addition, an analysis of the stability of the difference equations is presented 
and discussed, The numerical results are compared with experimental data 
for a 10-deg cone at 12-deg angle of attack, and a 5, 6-deg cone at &=deg angle 


of attack. The agreement is very goou, 


A description of the associated computer program is contained in the 


appendices. 
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Re 


u,v, Ww 


x,Y,® 


SYMBOJS 


static enthalpy (h = n/n) 

finite -difference grid points in x,y, directions, respectively 
number of mesh points in y and @ directions, respectively 
length used to nondimensionalize variables 

free stream Mach number 

local streamwise Mach number (M, = uM, /Nh) 

: : wy me 
dimensionless pressure (p = p/P VG) 
freestream Prandtl number (assumed constant) 
distance from a point in the flow to the axis of symmetry of 


the cone (r = x sin@ + y cos @) _ 
eV y’\ 


Ww UW 
~ 


free stream Reynolds number 7 
Boo 


Sutherland constant 


dimensionless velocity components in x,y, ® directions, 
res i =u/V_,vev/V_, w=W/V 
espectively (u i a / S wiv 

free stream velocity (V_ = V_/V¥_ = 1) 


ag ao co 


coordinates along the cone, normal to the cone, around the 
. 7 . owen tae eas 
cone, respectively (xX = Xi¥, y = Tit) 


angle of attack 

ratio of specific heats (assumed constant) 
mesh spacing in x,y,® directions, respectively 
transformed normal coordinate (n = y/£) 

cone half angle 


viscosity [bh = Vh (1+S)/(1+S/h)] 
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SYMBOLS (Continued) 


distance from shock to cone surface 


? 
density (p = YM, p/h) 


denotes free stream conditions 


denotes conditions at the cone 


denotes dimensional quantity 


denotes iteration number 
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SECTION I 


INTRODUCTION 


In recent years there has been a tremendous incre..e in numerical 
solutions for three-dimensional flow problems. This increase has been due 
to the rapid growth in the storage capacity and speed of computers, The 
primary effort in computing three-dimensional flows has been in using 
explicit methods, Explicit methods, although relatively easy to program, 
consume excessive amounts of computer time due to stability restrictions on 


step sizes. Even a DuFort-Frankel (Ref. 1) or Crocco (Rei. 2) scheme uses 


a considerable amount of time when the mesh spacing is small, as is neces- 
sary for accuracy with high speed, high Reynolds number laminar flow prab- 
lems. Therefore, for many flow problems it is desirable to use an implicit 
technique to solve the governing partial differential equations. Implicit 
methods have the advantage of being stable, consistent, and accurate for 
reasonable stepsizes. The major drawback is the size and complexity of the 


computer program which must be written and the storage requirements due 


be ete een me GOI jaa nla lla etic 


to the necessity of solving large systems of equations. Alternating direction 
implicit methods (Ref, 3), although reducing the size of the system of equa- 


tions which must be soived for three-dimensional problems, double the 


ee toe eared 


complexity of the computer code which must be developed. 


in this report, a method of solving the large system of algebraic equa- 
tions which result from the implicit ditferencing of three-dimensional flow 
equations is developed. For atypical problem, the implicit differencing may 
result in a system of 6000 ox more algebraic equations. A method of solving 
these equations which does not require excessiv2 coniputer storage and that 
yields accurate results is presented. The method is similar to the "predictor 


corrector" rultiple iteration technique described by Rubinand Lin (Ref 4), 
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The numerical approach begins with en implicit differencing of the 


systern of nonlinear partial differential equations, The nunalinear algebraic 


equations resulting from this differencing are first linearized, and the 
resulting linear algebraic equations are then solved using a Gauss-Seidel 
(Ref, 5) iteration method, The details of the method are given in Section II 
fer a sirmple model equation. Because of the necessity of iterating, which 
results from the numerical avproach, the convergence of the iteration 
snould be considered. This question is also analyzed in Lection Il for the 


model equation. 


In Section Il, the numerical technique which is developed is applied to 
the solution of an approximate system of three-dimensional equations which 
has been developed to predict the flow fields around cones at high angles of 
attack. -This approximate system of viscous equations has becn derived 
from the steady Navier-Stokes equations by assuming the gradients of the 
shear stress in the streamwise direction are mich sinalier than the gradients 
in the normal and circumferential directions (Ref. 6). The reswiting equa- 


tious are similar ta those developed by Lin and Rubin (Ref. 7) to solve the 


sharp tip, Jow Feynolds number region ier a cone at angle of attack. The 
resulting system of equations is first order in x and second order in y and@. 


The convergence and stability of the system cf equations are discussed. 


Solutions to the system of equations are presented for two cases in 
Section IV. Yhe first, a 10-deg half angle cone at 12-deg angle of attack and 
a freestream Mach npamver of 8; and the second, a 5, 6-deg half angle cone 
at 8-deg angle cf attack and a treestream Mach number of 14,2, The condi- 
tions for these cases correspond closcly tu experimenta: data obtained by 
Tracy (Ref. 8} and Stetson and Ojdara (Ref, 9). Comparisons ef the nurmeri- 


cal resuits with the exnerimental data are shown. 
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SECTION II 


NUMERICAL TECHNIQUE 


In this section, the numerical technigue is developed and analyzed. 
To illustrate the approach, the followiny three-dimensional partial differen- 


tial equation is considered; 


2 2 
du du du au a-u _ 
oe ty Obs ae ee c, d20 (1) 


This equation is representative of the viscous flow equation. For three- 


dimensional steady tiow, a = v/u, b = w/ur, c = B/Repu, d = u/Repur“. 


j 


The following finite difference approximation formulas are used: 
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where Wt is the value of u at the grid point ik £ é 


a 


In adattion, in the more general case, a cross derivative term appears. 
For completeness the difference formula is defined, 


_a*u 


ane 7 Ha ea 


- f 
: (O41, kel, 2-1 Ws ied, £1) /4nde (3) 


The scheme for solving the differential equation is then completely 
implicit. To obtain Bid, xf (solution known at j) a linear system of equa- 


™W 


tions of order K * L must be solved, 


The method proposed to solve this systern of linear equations is the 
line Gauss-Seidel iteration method mentioned in Fox (Ref. 5) and Isaacson 
and Keller (Ref. 10). To be specific, suppose Weyl, x,t is a guess to the 
solution of the difference equations where n denotes the iteration number. 
wen the correction Ske which must be added to ra kf to give the 
solution satisfies, after rearranging 


Lets 1,22, 24\5 . tela 
7 2Bn” Ae j41, k=2,f as 2 2) “j41,k,£ *\ Zan 27541, k+1,£ 


a A® An / 
ju" jut due ae n ; 9% ni Se, 
* "Ox oN Oe 1 ef 2 ~ 200 \341,k, f+1 ~ 9541,k, 2-1 
an ae = J 27K,L J 2%, : 
pe (3 41,%, fel i THe (4) 


~——~ 


If the underlined terms on the right hand side of Eq. (4) are ignored and the 


resulting equations are solved in the order L- 1,2,...L, usiag the boundary 


conditions at f= 1, then the approxiniate eres denoted by U. j4l, kf should 
ntl 

be close to u, ie], k; f° Taking us 441, kf > wis, kf +: erage ae anew guess 

to the selution of Eq. (1) the wees is repeated to obtain us 2 and so0 on 


until convergence is achieved, 
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The above method has the advantage that L systems of order K must be 


‘ solved instead of one cf order L * K, This saves time and storage. 


A seeming disadvantage is that iteraticn is required, That is, solving 
L systems of order K just produces a "guess" to the solution of Eq. (4). To 


obtain a better "guess" the L systems have to be solved again, etc, 


However for nonlinear problems this is not a disadvantage. To solve a 
nonlinear system of equations, some form of linearization must be done 
(e.g., Newton-Raphson method) and then iteration is done to obtain an accu~ 
rate solution. The line Gauss-Seidel method may be used to solve the linear 
system, and instead of iterating to convergence the first iterate is taken as 
the next iterate in the nonlinear sequence of iterates, Experience has shown 
that the convergence of the nonlinear iterates is not severely hindered by 


not solving exactly for the iterate. 


Once the iterations have converged then the method is completely im- 
Ppiicit and so the single linear Eq. (1) is stable anc consistent. The primary 
question to be answered is whether or not the iterates converge. To con- 


sidex the cinvergence question write the difference equation [Eq. (4)] as 


yt 5 Ae grt} yt , PAX n ~yrti 
“uel, kA” ZBq \[j41, ke1,£ ~ “j4t, k-1,0) 7 208 \ je. k, £41 7 9441, 6, £-1 


_ esd yotl cs 2yntl 4 gat 
"ane \itl, kt jtl, kf 7 “j41,k-1,£ 
, / 
GgAx n 5 ntl ntl s 
cae \ tleke del Yee (ieee bay. tied @) 
© j 


The above equation is a diflerence equation with difference index n. To 


determine under what conditions the solution converges as n-+o the Fourier 


series method as presented in Richtmyer and Morton (Ref. ae be used, 
n+l i(m ,£4% + m,kAn) 
1 2 
for Lk ande 


, n —— ; wa eae 
for er ae The term Bod is ignered since it is independent of n, and 


That is, substitute very Ses m2k4t)) 


the resulting equation is solved for A. The iterates will converge if |xki < 1. 


For Eq. (5) the amplification factor \ is 


Fa (ss si us) (cosm Ae + isinm 
Ae 


o*) 
N= BAx  dax cax ,, edax aAx BAx  dAx 
bias Gox ., a , ee 8: ae oe \ soe 
. 1- (sx +e ) coam ae ae (costa,An 1) + hee + re isinm,4n + (a3 + 2 | isinm Se 
(6) 
i : ; : A® 
Atter some manipulation it can be shown that IAI < 1 ifax < Tbr’ 


If the values of u at £-1 were evaluated at the a iterate instead of the 
anh: this would be the line Jacobi elimination method (Ref. 10) and the con- 
vergence criterion would be the same, It is line Jacobi elimination that 
Rubin and Lin (Ref. 4) studied; however, instead of considering convergence 
they looked at wat would happen if just one or iwo iterations were carried out. 
They found that the equaticns were not quite consistent and for the method to 


be stable as a marching scheme in x there was 


a restriction uudx depending 
on the number of iteratiot.s performed. For one iteration, the stability 
restriction was the same as the above derived convergence restriction. 
(Rubin and Lin considered *he case where c = 0,) The question may be a.ked 
why it is necessary to iterete to convergence (which may require four or five 
iterations) instead of iterating only once or twice. Iterating to convergence 
produces consistency, and in nonlinear equations it is necessary to iterate 
several timcs to obtain an accurate and stable solution to the nonlinear dif- 


ference equations. 


Tie rates of convergence for the line Jacobi method and the line Gauss - 
Seidel method have been studied for elliptic problems (Ref. 10) and it has been 
found that the line Gauss-Seidel method converges twice as fast as the line 


Jacobi method, The results for the above simple parabolic case are analo- 


gous, 


Ir a DuFort-Frankel scheme or a Crocco scheme or some cther modi- 


fied explicit formula (modified to remove the diffusive stability requirement) 


saith lamdhins 


+ ll all 


4 
is 


}! 


sated 
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is used, there are two convective stability requirements. The above scheme 
eliminates one convective Ax restriction. For problems which permit very 
unequal meshes in the two directions, as for many flow problerias, this may 


permit much greater step sizes, 


In the next section, the differencing described above is applied to a 
complicated system of three-dimensional viscous flow equations which have 


been developed to solve for the flow field around a cone at angle of attack. 


St 
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SECTION IIL 


APPLICATION TO CONE AT ANGLE OF ATTACK 


A. GOVERNING EQUATIONS 


The numerical technique which was developed in the previous section 
will be applied to a complicated system of three-dimensional viscous flow 
equations which have been derived to predict the flow around a cone at angle 
of attack, The system of equations has been derived from the steady Navier- 
Stokes equations by assuming the gradients of the shear stress in the stream- 
wise direction are much smaller than the gradients in the normal and cir- 
cumferential directions (Ref. 6), Yhe coordinete system used in the develup- 


ment of the equations is illustrated in Figure l. 
The resulting nondimensionai equations are listed below: 


Continuity equation 


dgpur , Opvr. apw _ g 
ax ' ay * ao * (7) 
x-rnomentum equation 
opu r aouvr aown De oe 35 
5x aT a + “3a - pw sinO® +r Dx 
> =e (of) + 2(2es).4 i 
“ Re | dy ey : re 35 (H ity ay cos 6 
(8) 
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= 
y-mMomentum equation : 
{ 
opuvr apver apvw ep : 
9x ey rs - pw cos® + T by 
ee eos re pe gan (te) : 
= Re | 3 ay \" by * 2 a6 \“ ax 
i 
1 2 (,a@)_ 2 2 fu, ow : 
aes ) Far le ae) (9) 4 
?.-momentum equation 
opuwr opywr a pee ap 
ae Rae “ay “spt puw sin®@ + pvw cos®@ + aD 
ro (ev) Be (ue), 2 (uae), 4 8 (aw | 
Re | ay C ag/ ~ 3x ag \P ay)? ay a ak gree (u 28 )| 
(10) 
Energy equation 


(11) 


where the perfect gas equation of state is used to relate the density to the 
pressure and enthalpy 


p= YM* P (12) 


and Sutherland's law is used to relate the viscosity tc the enthalpy 


148 
p = vh TtS/h (13) 


A censtant Prandtl number and specific heat will also be assumed. 


The above equations are similar to those used by Lin and Rubin (Ref. 7) 
except the terms associated with r—0O have been dropped. These terms are 
important only near the tip at low Reynolds number. We will be interested 
in solving the higher Reynolds number cases downstream of the tip region, 


The following boundary conditions at the cone surface are used, 


u=vewoe20 

h. = specified constant (14) 
jae, 2 (4 av, 1 p atw) 

\ay/ Re 3 oy r 3 ayag / 


The last equation has been obtained from the v-momentum equation using the 


condition (av/dy) = 0 which is required in order that the continuity equation 
is satisfied at the wall. 


The following difference formulas are used in the (ep/ay), equation: 


| and are obtained from Taylor series expansions in the normal direction of 
aw 


OP lh. -2 k=2° 
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The Rankine-Hugoniot jump conditions are applied at the shock boundary. 


In the body-oriented coordinate system (Figure 1) they are 


—'s 


Conservation of mass equation 


ag 1 9f 
= ) a = - = r . — , — —— = 
(Ua 7 PKG By (Mase PVR) 4 (Wis Pr Wy) = Bo 0 (15) 
Conservation of normal momentum equation 
(o,, #-v, rw 2 aby (a Poca ee LOY 
oo Ox cor o@¢ Ep seni if K ax K Kr ae (16) 
2 2 oOo ° “KC 2 : 2 
ag 1 ag ag 1 0g§ 
(22) sie(4 ae ) =.) sae: =.) 


[ st eecel i = 
1 38 @ ' 1 96 86 _ 
(uo - ee + (= 3} | 4 We. - Vic) ae (We a Wy) = 3 ae = 0 
1 a€ 
Conservation of energy equation 
2 2 
(Y-1)M (Y-1)M 
2 v2 2 2 2 


The subscript K denotes the value of the variable just inside the shock. 


order to uniquely determine the six unknowns €, 


In 
bad 4 Vue Wie? Pre hy, the 
above five equations must be augmented with a sixth equation, A one-sided 
differencing of the continuity equation provides the sixth equation, Full 


justification and discussion of the above equations and boundary conditions 
are presented in Ref, 6), 
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Since the fluid flow is symmetric about the plane @ = 0 and @= 1, the 


equations used will be solved for 0 s<®st where the symmetry conditions 


ws Jps Age 7 Q 
aD (u,v, p,h, €) Sead 0; woe = 0 (19) 


are used at@=Oand@=t. (Note: ® = 0 is the windward side. ) 


The shock distance is to be solv-d for from the Rankine -Hugoniot jump 
conditions, Ina rectangular y-@ grid the shock may not fall on a mesh point 
so that mesh points would have to be moved or added to accommodate the 
shock, Thus, the transformation 9 = y/é (x,@) is made. The resulting equa- 
tions are then solved for OS yn sl, 0£ @<s 1, where yn = 9 corresponds to the 
cone and n= 1 corresponds to the shock, The shock distance € appears in 
all the equations and in order to keep the matrix of coefficients obtained from 


the difference form of the equation in block tri-diagonal form, a sixth equation 
Ca (20) 
on 


is differenced. Thus the problem to be solved consists of six differential 
equations, six boundary conditions at each of the positions yn = 0 and n = 1 and 
two symmetry conditions at @ = 0 and @=w [Eqs. (7) through (20)] in the six 
unknowns u,v,Ww, p,h, &. 


If initial conditions were known, a marching scheme in x could be used 
to solve the equations. An explicit method would not work since, as shown 
by Baum and Denison [Eq. (12)] for the axisymmetric problem, it is not pos- 
sible to solve for du/8x, av/ax, AwW/Ox, Bp/ex, dh/dx, at M,, =l. Evenif 
this difficulty could be overcome the diffusive stability requirement would be 
too strict, Since the gradients in the normal direction are much larger than 
in the circumferential direction, the normal mesh will be mech finer than the 
circumferential mesh, Thus a method that is implicit in the normal direc- 


tion should be more efficient than a modified explicit differencing such as 


DuFort-Frankel or Crocco, Accurate solutions at each x-station are desired 


and since the equations are very nonlinear they require iteration for accuracy. 


The method proposed in the previous section is most appropriate fog this 
problem. 


An alternating direction implicit technique was tried, with iteration to 


handle the nonlinearities. It was found that near u = 0 and M.. = 1, the equa- 
tions for the implicit in @ step were ill-conditioned and meaning[ul sclutions 


could not be obtained. The M,. = 1 difficulty tvas overcome by evaluating 


the 6p/dx term backwards in x in the x-mormentum and energy equations. 


However, the u = 0 difficulty reimained, 


The implicit difference equations |Eqs. (2) and (3)) are substituted into 
the partial differential equations and « system of nonlinear algebraic equa- 
tions result. There are many ways to linearize such a system. Since con- 
vergence is guaranteed provided the initial guess is close enough, and because 
the convergence is quudratic, the Newton-Raphson method is used iv solve 
these equations, That is, the nonlinear terms are expanded in a Taylor 
series and terms higher than tirst order are dropped. It is known that this 
iteration procedure converges provided the initial guess is close enough to 
the solution. It was found that linearly extrapolating the sclution at the pre- 
vious two x stations gives a satisfactory initial guess, To see what lineari- 
zation does to various terme let f represent the increment to be added to the 


known iterate and the superscript n dencte that iterate. A few sarmple expres - 


rn ni n n 
gear _ (dpor (Pry), (ur) /8p>, dP \ , = (Pu) 
ax ( ox ) oa ax | Ax (225+ Shh) 2 4X 

a rn mn 
dew _ (apw\ 3p ws) a n Op - map 
fe) go + be (Wap tw BER) 
. vs ? n Pax 
pox. pay) +R(Sh 2 4 prey 
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The r term above is r= £ neos®, since Koay ain@ is Known. The linear sys- 


tein of equations obtained are solved for u,v,W,p.6, and € using the line 


Gauss-Seidel methed described in the previous sestion, However, cnly one 


iteration of the line Gauss-Seidel method is performed. ‘The approximate 


solutions so obtained wre then used to obtain the next guess to tne nonlinear 
system, . 


For this flow preblem the L systems of equations that must be solved 
to obtain one Gauss-Svidel iterate are of order 6 * K since there are six 
variables involved. The matrix of coefficients is of block tridiagonal form, 
An efficient methad for solving such systems is presented in Isaacson and 


Keller (Ref. 10) and was used by Rubin and Lin (Ref. 4) and also in the 


present analysis. 
8. CONVERGENCE AND STABILITY 


The questions of convergence and stability niust be investigated for the 
system of flow equations. To consider convergence, the equations are simi- 
pliitied, birst the equations are written ina different form by expanding the 
derivative expressi: 43 and subtracting the continuity equation from the 
momentum equation and the energy equation, 


that 2/Pr = 4/3 pe=p. 


For simplicity it is assumed 
The viscosity is assumed to be constant and the 


shock distance is assumed to be known. Since the iteration is primarily for 


the ® derivative terms, then all terms not involving derivatives with respect 


tu @ exclusively are ignored. The following equations are left. 


gu, 1 ap, w dup * au 

8x Pu ax ur do Ke ag? et 

ov w Ov x ave 

a 1-——- g--Re SS = 0 (21a) 
x ur ¢® aac ‘ 

ow WwW aw x 3 Ww h ap 7 

ax . ur do R : sie 
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ae tan eee 


ce ant Ra tue 


conte UD nets EE 
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ox u h 


5x x ur eo ur 86 Ge eae 
422b) 
ent MOS 29 a P 
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| 1 uh 
Re = er .8}—=— eee 
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To further simplify the analysis, assume that all coefficients of the deriva- 
tive terms are constant, and that the mesh spacings are constantAx, Ay, and 


4@ in the x, v, and @directions, respectively. After differencing and come 


rearranging, the line Gauss-Seidel iteration method reduces the difference 


equation for wat ¢ach mesh point to 
n+l ok n+l pti Ax n 7 yn 
yal RL * bu Pjet, kL tur 2ae y41,k, b+) 7 4a, 6, f-1 


* Qx fon _ ayntl n+l _ 1 
Re Ape Gar Ms ee : “jt k £1) = Bet * PEP; xf 


(22) 


with similar expressions for the ether equations, where the superscript n 
cenotes the known iterate, the superscript n+l denotes the next (unknown) 


iterate, and the right hand side of the equations are at j and so are indepen- 
. cent of the iterate. - 


The convergeuce of the solution uf the system of difisrence equations 
represented vy Eq. (22) cam be determined using the Fourier series method 
as presented in Re’. 11, The solutions will converge if the eigenvalues of 


the associated amplification matrix are less than one in absolute value. To 


oe ye 


Scie Raith ai gee ee et had sek ot alll lt) sciatic dah 


: im£, irng 
obtain the eigenvalues substitute ue mbae and du eltnbat for wu? 


nel ee : : : ; i+, kL ae 
aL, cf! and similarly for v,w,p,h; this gives five linear equation in 

ur Yor Wo: Po and ho. and the eigenvalues of the amplification matrix are 
those values of \ for which the determinant of the coefficient matrix is zero, 

For the above system it has been found in regions of the flow where 

M, is near one or less than one, that some of the eigenvalues of the amplifi- 
cation matrix are greater than one, in absolute value, independent of Ax and 
Ab. However, if the ap/ax term in the x-momentum and energy equation are 


differenced backwards in x, i.e, 


op = Ps, k£ ° Pi-1,kb 
ox Ax 


(23) 
as suggested by Ohrenberger and Baum (Ref, 13), or set equal to zero as 
suggested by Rubin and Lin (Ref, 14}, then convergence criteria can be ob- 
tained, Doing either of the above modifications causes the difference expres - 
sions representing the underlined terms in Eq, (21) to be independent of the 
iteration index n. Thus, the first and second difference equations are 
uncoupled from the rest of the system and are similar to the simple equation 
in Section II, Therefore tai < lif 


Ax < |[=lae (24) 


The remaining three equations are studied using the Fourier series 
method, If the determinant of the associated three by three matrix of coef- 
ficients is set equal to zero, it can be found that a third value of ) is the 
game as the first two above, and so Eq, (24) must be satisfied. The remain- 
ing two vaJues satisfy a complicated quadratic equation. To simplify the 
analysis the first and second derivative terms are handled separately. For 


second derivative terms only, it can be easily shown that both roots of the 


iy ine eat oe eae 


quadratic equation are less than one in absolute value, For first derivative 


terms only, it can be shown after considerable effort that the two roots of 


the quadratic equation are less than one in absolute value if and only if 


RR <r EE (28) 


Iwt (1 + 1) 4 Vl iy we + =a 


The results as represented by Eqs, (24), (25) apply to the simplified 
linear system Eq. (21). It has been found numerically that the restriction 
for the actual nonlinear system of equations [Eqs. (7) through (20)] is quali- 


tatively like Eqs. (24), (25). Quantitatively the restriction is similar to 
Eq. (24). 


Once the solution can be obtained at Xia! the question of whether the 


scheme is stable for marching in the x direciiun must be answered. To 
analyze the stability, the following system is considered. ; 
3 

au 1 3p * au 2a9u : 
dx * pu ax Zz 2) * ; 
co) ay i 

eae | 

ene" (Ep 422 Hs). 0 : 
de ay 1 
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The above equations are differenced implicitly and the resulting egqua- 


tions are studied using the Fourier series technique. 
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The equations are 
stable if the eigenvalues of the amplification matrix are less than or equal to 
one in absolute value. 


ci, Anabel 


The second and third equations are uncoupled from 
the system. 
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It can easily be shown that they are unconditionally stable, 


For the remaining three equations one of the eigenvalues is A = 1, and 
another is 


1] 
Be ey ae 


epee One (27) 
* rAx : 
1~2 Re |(cosm,Ae-1) as + (coem,A-1) — 4 
1 me rd ‘ 
AY 4 


which is less than or equal to one in absolute value. 


Since the first four eigenvalues are less than or equa! to one in absolute i 


the stability of the system of Eq. (26) depends on the magnitude of the 


fifth and finalone. Three different cases are considered depending on how 


the ap/dx term in the x-momentum and energy equations is differenced, In 


value, 


Sores 


we 


the first case, if dp/dx is set to zero as suggested by Rubin and Lin (Ref. 14), 
then the [fifth eigenvalue is the same as Fq. (27). Thus, the difference equa- 
tions are unconditionally stable as a marching scheme in x, 


The second case corresponds to evaluating 9p/8x backwards in x 
[Eq. (23)], as suggested by Ohrenberger and Baum (Ref. 13). If M,, 


21 then 
it can be shown that the magnitude of the e ali 


igenvaliue is less than or ‘equal to 
one, However, if M,< 1, which occurs near the cone due to the boundary 
condition u = 0, the following restriction must be satisfied to ensure that the 


absolute value of the eigenvalue is less than or equal to one. 


Ax 2 (ae : ) —- | 
M,. 2YRe (1-cosm ,A®) 4+; + (1-copm,Ay) +5 
Ag Ay” 


(28) 
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The third case corresponds to taking 9p/ax implicitly. For this case 
it is found that Ax must be twice as big as for the previous case where 9p/dx 
was evaluated explicitly. This implies that if a different metLod were used to 
solve the implicit equations (the proposed Newton-Gauss~-Seidel method does 


not converge when dp/dx is differenced implicitly as mentioned previously) a 


numerical solution to the fluid flow equations could ve obtained ifAx were : 
chosen sufficiently large. A lower bound restriction on the marcning step- j 
size has also been found for certain stiff ordinary differential equat:ons by 3 
Curtiss and Hirschfelder (Ref. 15) in order to suppress so called de:arture 
solutions. The case here is analogous to the ordinary differential eguation 3 


case. The departure ig characterized by the leeward surface pressure 


oscillating or rapidly increasing, and has been observed by Baum and Denison 
(Ref. 12), Rubin and Lin (Ref. 14), and Tyson (Ref. 16). 


Tyson experimen- 
tally found that a large stepsize was necessary to suppress the departure 
‘solutions. ; 7 
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SECTION IV 


NUMERICAL RESULTS 


To demonstrate the validity of the technique, solutions have been com- 
pared with experimental data obtained by Tracy (Ref. 8) on a sharp 10 deg 


hal’ angle cone at an angle of attack of 12 deg. The parameters used in the 


These correspond closely to the experi- 
Initial conditions, which would normally come from a solntion 


calculation are given in Table I, 


mental data. 


io the nose region, are needed before an exact comparison can be obtained, 


Since the nose solution was not available, the following technique was used to 


generate the required initial conditione. Starting at zero, the angle of attack 


was slowly increased while marching along the cone until 12 deg was reached 
atanx=xX.. 
fe) 


The calculations are then continued at a constant 12 deg angle 


of attack, and the solution is allowed to relax to the desired sharp cone 


results. Because of the method used to generate the initial conditions, the 


calculations are not expected to agree with the data at the same x station. 


Symbol and Value 


—er 


i 

1 

Yowever, the calculation should relax to results which are similar to the 4 

: data (except the difference in local Reynolds number) as the solution continues 4 
downstream, | 

ae Table I. Parameters for Tracy's Case | 3 

r Paraineter i; 


eone half angle 


2 of attack 


Greestceam Reynolds number 
freeatrearn Prandtl number 


{ 
| reiio of specific heats 
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Figures 2, 3 show the experimental surface pressure and heat transfer 
arcund the cone at x = 0,33 it and the calculated results at x = x/x, = 8.5, 
25, and 50. The calculated heat transfer results illustrate that at x = 8,5 the 
effect of tlh: initial conditions at x = X, have not yet relaxed. At x = 25 and 
50 the caiculations appear to be approaching a relaxed result and the agree- 


ment with the heat transfer data is very good, 


In Figure 4, a comparison between the measured and calculated bow 
shock and viscous layer thickness around the cone is given. The data were 
taken at x = 0,286 ft, the calculated results at x = 50 are shown. Finally 
in Figure 5, the calculated velocity vectors on the leeward side projected 
normal to the streamwise direction are shown. The separated flow region on 


the leeward side is clearly shown by this figure. 


The above results were obtained by evaluating dp/dx explicitly. Runs 
were also made with gp/ax = 0 and the results differed very slightly from 
those presented. 


A comparison of actual convergence and stability restrictions with the 
analytical restrictions for this case is presented in Figure 6. From Eq. (24) 
with A® = 10 deg and taking 1 as a lower bound for tu/w| we see that 


ax 
a <0. 03 (29) 


is yeCeseary for convergence. Equation (28) implies that the biggest restric- 
tion occurs for small u, that is near to the cone. To simplify the expres - 
sion, suppose that cos mj A® and cos m,4y are zero, This is not strictly 
valid; however, we are looking for an approximate answer. Making this 
assumption may not give us precise quantitative results but hopefully quali- 
tative results will be obtained. For the cases run, Sy << r4é The nearest 


point to the cone that is solved for corresponds to y = Ay. The functions p 


and » evaluated at Ay are approximately constant as functions of x, Also it 
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Figure 3, Circumferential Surface Pressure Distribution 


for Tracy's Case, @~- 12 deg 
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Figure 4, Circumferential Heat Transfer Distribution 
for Tracy's Case, @ = ]2 deg 
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has been observed that the ratio Ay/u where u is evaluated at y = Ay is almost 


constant as a function of x. Thus Eq. (28) is approximately 


ax . [Se ax | Ay. (30) 
x pan x 

where the expression in brackets in Eq. (30) is approximately constant, In 
addition, since the transformation y =1€ is made, Ay =4%&; and since € 
increases approximately linearly with x, then Ay/x remains constant as x 
increases. Equation (30) then implies that for an initially chosen Ay, Ax 
must increase linearly as x increases, Therefore if Ay is chosen 80 that 
Eqs. (30), (29) are satisfied initially, they will be satisfied for all x p~ 


vided Ax increases proportionally to x. 


Tracy's case was run with Ax/x = 0.03 and a An spacing such that 
initially Ay/x = 0.0006, The expression in brackets in Eq, (30) is equal to 
lé so that Eq. (30) is satisfied. With S&= 10 several different values for 
Ax and Ay were used to determine how accurate the derived inequalities are. 
The results are shown in Figure 6. The shaded region is the predicted area 
where convergent stable solutions should be obtained. The actual region of 
convergent siable solutions is somewhat different, However, the qualitative 
results are correct, It was found that increasing Ax slowed down the con- 
vergence of the iterative procedure to obtain the solution at an x station until 
finally it did not converge. Decreasing Ax sped up the convergence but 
produced solutions that departed as a function of x, Increasing Ay led to 
departure solution and decreasing Ay led to stable solutions. For a few cases 
the results for the surface pressure on the leeward side (the most critical 


for departure) at a = 12 deg are shown in Figure 7. 


A second case has been run to compare with experimental data obtained 
by Stetson and Ojdana (Ref. 9) on a sharp 5.6 deg half angle cone at an angle 
of attack of 8 deg. The parameters used in the calculation are given in 
Table Il. Figure 8 shows the wall pressure distribution on the leeward side 


atx = 17. Stetson's case was run with ép/ax = 0. Analysis of the restrictions 
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“ on Ax for @p/3x evaluated explicitly as was done for Tracy's case led to Ax 

. spacings that were much too small in the sense that too many y mesh points 
were required and the storage capability of the computer was exceeded. 
Numerical experimentation verified that the solution could not be obtained 


for any permissable mesh spacing when ap/dx was evaluated explicitly. 


= Table Il, Parameters for Stetson's Case 
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cone half angle @ = 5,6 deg 

angle of attack a = 8 deg 

freestream Mach number M., = 14,2 ; 

freestrear. Reynolds number Re = 0,83 x10 /ft 

freestream Prandtl number Pr = 0.75 A 
ratio of specific heats v = 1,4 . 
Sutherland constant S) = 4 . 
freestream dimensionless pressure Pai = 0, 00354 


static enthalpy at the cone ho = 32,398 
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SECTION V 


COMPUTER PROGRAM 


A computer program has been developed for the CDC/760( computer 
that solves the above equations, Provision has been made for v. .'1ble grid 
sizes in both the y and @ directions and this feature has been used heavily. 
The iteration logic has been structured so that when the solution along a® 
line converges to the desired number of figures, that line is dropped from 
the iteration loop. This saves considerable time since a few lines require 
as many as seven iterations to converge to five figures while most @ lines 


require only two or three iterations, 


For 50 points in the y direction and 23 points in the @ direction the 
program requires 35,000 words of storage, It takes 30 sec to obtain the 
solution accurate to five figures at one x station. Eleven percent of thi 
is spent evaluating all the derivative expressions, 33 percent is spent com- 
puting the Jacobian coefficient matrix, and 56 percent is spent solving the 
linear equations. For Tracy's case it took about 120 steps to go from 
x = To RS 50, 


For convenience in programming the derivatives were defined as 
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The right hand side of Eq. (4) became 
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where the underlined terms in Eq. (4) have been ignored. It is seen that 


it is necessary to store only two planes of the solution, one at j and one for 


: the current iterate at j+l. The n+l iterate is stored on top of the n iterate 
as it is computed, . 
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SECTION VI 


. DISCUSSION AND CONCLUSIONS 


A method for solving the implicit difference equations describing the 
three-dimensional flow around a cone at angle of attack has been described 
and analyzed, These equations which are derived and discussed in Ref. 6 
constitute a system of three-dimensional nonlinear parabolic equations. The 
technique for solving these equations has been shown to be accurate and 


efficient in both running time and computer stcrage. 


The numerical method is not restricted to steady flow problems but 


could easily be applied to two-dimensional time dependent calculations, 


An analysis of the numerical aspect of departure solutions has been 
presented, Methods that have been proposed by other authors to suppress 
the departure solutions have been verified and in some cases qualified, In 
addition, it has been shown that departure solutions can be suppressed even 


if the streamwise pressure derivative is included, if the steps ze is large 
enough, 


Results have been compared with experimental data and the agreement 
is very good. 
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APPENDIX A 


PROGRAM INPUT INSTRUCTIONS 


The input is divided into two sections, namelist input and formatted 


Parameters describing the problem are read using NAMELIST/INPUT/, 


GAMMA = Y 
MINF = M 
@ 

THETAC = @ 

REINF = Re 

PRINF = Pr 

ALFA = @ 

PINE pe 

SPROP = § 

NJ = number of x-stations 

NK = number of y-stations at which initial conditions are 
specified 

NL = number of $-stations at whic! initial conditions are 
specified 

MOD = .TRUE. or . FALSE. 
depending on whether the input rnesh distribution is to be 
modified oz not. Default is . FALSE, 

ITAPE = 0, no output on TAPE2, this is the default 


= N, output solution on TAPE2 every Nth x-step 
The rest of the input is read using format control, 


IREAD (112) 


IREAD > 0 read initial profile from TAPE3, used for restart 
purposes, The first profile on TAPE 3 at an x station 


bigger than or equal to X(1) is selected as the initial 
profile. 


IREAD £ 0 read initial profile from cards 


pees te AE 48 


les ere Sere 


<> RAR EE 


CEE LGR en Sea ne ae et hia tel Reiterates da hae Bae Mapes Reet e i2 ff cae Re gn ad ete ERTS 


If IREAD s 0 then the cards with the solution are input 
FI(L), ZI(1, L) (2E12.5) 
repeat FI(L) is a ® station in degrees 
ve sme e Z1(1, L) is the shock distance at @ 
. ET(K), U(K, LL), V(K, L), W(K, L), PUK, L), H(K, L), Kel, NK(6E 12, §) 
F,T(K) is a distance from the body (=y not y/shock) 
U,V,W, P,H are values of u,v,w,p,h at ET(K), FI(L) 


X(J) J=l or J=1, NJ (6E 12.5) 
%(J) is an X-station along the body, If J=1 then the rest of the 
points are obtained using X(J!1) = c * X(J) where 


c is a constant > 1, X(1)is the position of the initial profile. 


The following cards describe v and h at the body: 


12) 
Lis the index of a ® station 
' J, VB(J,L), HB(J, L) (112, 212.5) 
J is the index of an X station 
VB(J, L) ic the value of v at the body at X(J), FI(L) vag 
HB(J, L) is the value of h at the body at X(J), FI(L) 4 
This card is repeated with N increasing from 1 to NJ, : 
tf any J stations are skipped then linear interpolation is | 
used to obtain VB and HB at the skipped stations 
There then follows a card with another L value and cards with J, 
VB, and HB. These groups are repeated with L increasing from 


lto NL. Ifany L stations are skipped then linear interpolation 


is used to obtain VB anc HB at the skipped stations, 


If MOD-. TRUE. more cards are needed. 

NEWK, NEWL (215) 
NEWK is a new value tor NK, and indicates that the input 4 
distribution is to be changed. If the n distribution is not to 
be changed then leave NEWK blank. 
NEWLas with NEWK but for NL and @. 

If NEWK > 0 
ETNEW(K) K=1l, NEWK (6D12. 4). 
The numbers must increase from 0. The input numbers 
are scaled by the program to go from 0, to 1, by dividing 
by ETNEW (NEWR) 

If NEWL > 0 
FINEW(L) L-1, NEWL (6D12, 4) 


These are the new @ stations in degrees, They must 


increase fram Q to 180 deg. 


ee ns 


APPENDIX B 


PROGRAM OUTPUT 


The nainelist input is printed out. Ifthe initial profiles were read 


fron cards then the card images are printed. 
The x-stations at which the solution will be obtained are printed, 


The input profiles are printed and if the mesh distribution was changed, 


then the new initial profiles are printed, 


The profiles are printed at each x-station as they are obtained. Pre- 
ceeding cach station printout is printed the iteration and convergence history 
for that x-station. The variables printed are ITER and INS(L), Lo], NL. 

If all @-rays have converged then ITI.R=-0, otherwise ITER:1]. Ifthe Lth 
@-ray has converped then INS(L) 0, otherwise INS(L) = 1. 


The parameters read in with namelist are written on TAPE4, followed 
by the profiles at each x-station. The data are written without format coen- 


trol and can be used te supply initial conditions. 


If (TAPE > 0 then the profiles are written on TAPE2 at each ITAPEth 


x-station, TVAPE2 can be set up as the punch file, 


APPENDIX C 


EXAMPLE PROBLEM 


The input necessary to run Tracy's case (sce Table I) is given in 
Appendix D. The sequence numbers on the initial condition profiles are not 


required, they are included on the sample deck to aid in ordering the cards if 
they get mixed up. 

The output from the program is presented in Appendix E. The solution 
has been printed at only a few mesh points. 


After running much further along the cone the solution profiles obtained 


would be similar to those shown in Figures C-1l, C-2, C-3, C-4, and C-5, 
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Figure C-1, Streamwise Velocity Profiles for 
Tracy's Case, a = 12 deg 
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Figure C-2, Normai Velocitv Profiles for 
Tracy's Case, a= 12 deg 
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Figure C-4, Pressure Profiles for Tracy's Case, 
a= 12 deg 
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Figure c-5, Enthalpy Profiles for Tracy's Case. 
a= 12 deg 
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APPENDIX D 


INPUT FOR TRACY'S CASE 
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